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The relaxed and unrelaxed formation energies of neutral antisites and interstitial defects in InP are calculated 
using ab initio density functional theory and simple cubic supercells of up to 512 atoms. The finite size errors 
in the formation energies of all the neutral defects arising from the supercell approximation are examined and 
corrected for using finite size scaling methods, which are shown to be a very promising approach to the problem. 
Elastic errors scale linearly, whilst the errors arising from charge multipole interactions between the defect and 
its images in the periodic boundary conditions have a linear plus a higher order term, for which a cubic provides 
the best fit. These latter errors are shown to be significant even for neutral defects. Instances are also presented 
where even the 512 atom supercell is not sufficiently converged. Instead, physically relevant results can be 
obtained only by finite size scaling the results of calculations in several supercells, up to and including the 512 
atom cell and in extreme cases possibly even including the 1000 atom supercell. 
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I. INTRODUCTION 

Over the past decade or so first principles density func- 
tional theory (DFT) 1 has become a powerful tool for study- 
ing the properties of defects in semiconductors. It is possi- 
ble to calculate formation energies, binding energies and mi- 
gration barriers, to predict local structure and, up to a certain 
point, defect levels and electrical activity. Problems and lim- 
itations remain, however, and one of the greatest is the very 
limited size of the systems for which calculations are feasible: 
10s or 100s of atoms, even when we wish to describe phys- 
ical problems involving 1,000s or 10,000s. This leaves the 
results heavily influenced by errors arising from the bound- 
ary conditions. These errors must therefore be carefully stud- 
ied so that their effects can be understood and accounted for 
when results are interpreted. There are two types of boundary 
conditions commonly used: open and periodic. Open bound- 
ary conditions are usually encountered in cluster calculations. 
The surface atoms are "terminated" with hydrogen to use up 
spare electrons, but are otherwise surrounded by empty space. 
Periodic boundary conditions (PBCs), meanwhile, are found 
in supercell calculations, in which a block of atoms is sur- 
rounded not by empty space but by an infinite array of copies 
of itself. Both approaches have strengths and weaknesses. We 
present here a detailed study of the problems arising from the 
use of the supercell approximation and we propose a method 
to overcome them. 

Finite size errors in supercell calculations come from two 
main sources. Elastic errors often arise because the supercell 
is simply not large enough to contain all of the local relax- 
ation around the defect, leaving calculated formation energies 
too high. In addition, the defect interacts with an infinite ar- 
ray of spurious images of itself seen in the PBCs, via both 
elastic and electrostatic interactions. The direct elastic inter- 
actions can easily be truncated by freezing all atoms lying on 



the surface of the cell at their ideal lattice positions. The elec- 
trostatic interactions, on the other hand, cannot be truncated 
or removed. They result in errors in the calculated forma- 
tion, binding and migration energies, errors which can be on 
the same order as the energies themselves. For practical su- 
percell sizes they need not even be negligible for neutral de- 
fects, since dipolar and quadrupolar interactions can remain 
significant. These latter can even result in errors in the cal- 
culated structures, since they favour certain symmetries and 
local relaxation modes over others. Hence indirect elastic er- 
rors cannot be avoided either. Finally, a third source of finite 
size errors is also present: the defect state wavefunctions can 
overlap with their images in the PBCs leading to a spurious 
dispersion of the defect levels which in turn can affect the for- 
mation energies, especially if the defect level is only partially 
filled. The errors related to this dispersion (or tunnelling) are 
expected to have only a fairly small and rather short ranged 
(exponential) effect. 

Recently, various correction schemes have been 
suggested 2 - 3 to compensate for at least the leading terms in 
the errors arising from the electrostatic interactions. They 
are usually based upon fits to quasi-classical models and/or 
multipole expansions of the electrostatic interactions. They 
have met with varying levels of success but, so far at least, are 
generally considered insufficiently reliable. There are more 
direct approaches, however. Probert and Payne 4 recently 
presented a detailed ab initio study of the neutral vacancy 
in Si, considering all aspects of convergence, from basis set 
and k-point sampling to size and symmetry of supercells. 
They demonstrated that the use of "large" supercells (200+ 
atoms) can be essential for obtaining the correct physical 
results. Meanwhile, we recently presented 5 a study of the 
neutral vacancies in InP. We demonstrated the advantages of 
not only using large supercells but also finite size scaling the 
results obtained. We evaluated both the relaxed and unrelaxed 
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formation energies of the phosphorus (Vp) and indium ( 
vacancies in simple cubic supercells of 8, 64, 216 and 
atoms. We then showed that the variation in the forma 
energy with supercell size, L, follows rather closely the fc 

E^{L)=E^ + a l L- 1 +a 3 L- 3 

where E^(L) is the formation energy in supercell "c" an 
and «3 are fitting parameters. is the finite size scaled 
mation energy corresponding to an infinitely large super 
Equation ^ has, in fact, the same form (L 1 plus L 3 ) as 
corrections proposed by Makov and Payne 2 . We will retui 
the issue of whether or not this is the correct form in sec 
IIVI For the vacancies we showed that the error bars on 
values obtained for E°f are usually rather small: O(0. 1 e\ 
less, depending also upon the level of k-point convergenc 
individual cells. Care must be taken, however, when SCaiini; 
relaxed formation energies of strongly Jahn-Teller 6 active de- 
fects, such as Vp. In such cases rather wider error bars are 
obtained if the symmetry of the relaxed structures varies with 
supercell size. Numerically, we also found that (for example) 
E°f for the unrelaxed Vi„ is actually ^0.2 eV above the value 
obtained in the 512 atom supercell, suggesting that there are 
cases for which even the 512 atom cell is not large enough to 
be converged, so that scaling becomes essential. 

In the present paper we extend the study to that of the other 
neutral native defects: the antisites and interstitials. It should 
be noted, however, that our primary purpose is not the study 
of the defects themselves but of the finite size errors which 
arise when calculating their formation energies. InP has the 
zinc -blend structure, with two antisites, Pi n and Inp, both with 
tetragonal (Td) symmetry when unrelaxed. For interstitials 
there are three high symmetry sites: two tetragonal, with four 
In or four P nearest neighbours (X^) or Xi ( p)) and a quasi- 
hexagonal site (Xi(hex)) with six nearest neighbours (three P 
and three In) and C3 V symmetry. Previous work for the iso- 
lated neutral cases has mostly been limited to the 64 atom su- 
percell. Here, formation energies were around 5-6 eV for the 
tetrahedral interstitials 7 and around 3 eV for the antisites^. Inp 
displayed some Jahn-Teller behaviour whilst Pi„ did not. We 
will examine how these results change with larger supercells 
and with scaling. In the next section we describe the method 
to be used. In section |III| we will describe the basic scaling 
results for the various neutral defects and will examine in sec- 
tions [W] and [V] the form of the scaling and when it does and 
does not work. In section lvTl we consider briefly other charge 
states of the defects. In Section IVm we discuss the origin of 
the surprising linear scaling we find in certain cases. In section 
IVIIII we estimate the size of the other non-finite size dependent 
errors for the purpose of comparison, before concluding in llXI 

II. METHOD AND K-POINT CONVERGENCE. 

We use planewave ab initio DFT— within the lo- 
cal density approximation (LDA) together with ultrasoft 
pseudopotentials 9 (US-PP) using the VASP code 10 . We re- 
cently presented 1 1 a study of the [Zni„-Vp] complex in InP us- 
ing the same technique and potentials. Unfortunately, physical 
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FIG. 1: Convergence of formation energy (N) with even NxNxN 
Monkhorst-Pack k-point grids, x : calculated values of formation en- 
ergies. Dotted line: "converged" value (calculated using 30x30x30 
k-point grid). Dot-dashed: running average over E9(N). Dashed 
line: running average weighted by the number of k-points in the irre- 
ducible Brillouin zone. Solid line: running average weighted by the 
total number of k-points in the full Brillouin zone. 

memory limitations in even the largest parallel computer facil- 
ities available to us mean that we must treat the 4d electrons 
of In as core, even though they are comparatively shallow. (A 
calculation for bulk InP in its FCC primitive cell places the 
In 4d states about 14.5 eV below the valence band edge at 
the r point when they are treated as valence - fairly deep but 
still close enough to potentially make some contributions to 
formation energies.) It would be preferable to treat them as 
valence, but then properly k-point converged calculations in 
the 512 atom cell - which our analysis requires - would be im- 
possible. The size of the resulting error will be examined in 
section lvTTll 

The optimized LDA lattice constant using these pseudopo- 
tentials is 5.827 A and the band gap is 0.667 eV, compared to 
5.869 A and 1.344 eV in experiment. As stated above, we 
will only use simple cubic supercells of 8, 64, 216 and 512 
atoms, since it is important to keep to a single supercell sym- 
metry since the scaling is different for different symmetries. 
No restrictions are placed upon the symmetry of relaxations, 
but we do not allow atoms located on the surface of the cell 
to relax. The exception is interstitials at the quasi-hexagonal 
sites in the 8 atom cell. Here, three of the nearest neighbours 
lie on the surface. Initial tests showed that the relaxation was 
not stable if all of these were allowed to relax at once so re- 
laxations were done in stages: firstly one group of neighbours 
was relaxed whilst the others were kept fixed, then vice verse. 

The key quantity is the formation energy 

E% = E$ (defect) - E$ (bulk) + £ mm (2) 

i 

where Ej (defect) and £j(bulk) are the total energy of the 
supercell with and without the defect, calculated using the 
same values of planewave cutoff, k-point grid, etc, to make 
use of the cancellation of errors. The defect is formed by 
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adding/removing n, atoms of chemical potential jU,-. We use 
Hp = 3.485 eV and pL\ n = 6.243 eVjii corresponding to stoicio- 
metric conditions. A planewave cutoff energy of 200 eV and a 
Monkhorst-Pack 4x4x4 k-point grid 13 was found sufficient 11 
for converged non-relaxed calculations in the 64 atom super- 
cell with errors around 0(0.01 eV). Hence grids of 8x8x8 in 
the 8 atom cell and 2x2x2 in the larger cells should suffice. 
However, in this study we examine specifically the errors aris- 
ing from the supercell approximation itself. K-point conver- 
gence is different in different supercells and we do not wish 
to include any significant errors due to this. Hence we need to 
ensure even higher levels of convergence during this study - 
much higher than is normally required or practical. (Examin- 
ing the scaled properties of the native InP defects themselves, 
over all relevant charge states, is left to future work, however, 
as it will be done with rather different levels of k-point con- 
vergence, pseudopotentials etc, once we have in the present 
work established the behaviour of the finite size errors.) We 
use k-point grids of up to and including 12x12x12, 8x8x8, 
4x4x4 and 4x4x4 in the 8, 64 216 and 512 atom cells respec- 
tively. (Exceptions are: Pi n - only 2x2x2 needed in the 512 
atom cell, and In;(i n ) and In;(p) - 6x6x6 used in the 216 atom 
cell to check convergence.) To improve convergence further 
we use weighted averages over £ d values calculated using a 
series of grids. Fig^shows the advantages clearly: the unre- 
laxed formation energies Ej° s (N) for Pi(i„) in the 8 atom cell 
are shown, calculated using NxNxN k-point grids and plotted 
against N. (This case has the most difficult k-point conver- 
gence in this paper.) To get errors safely below 0(0.005 eV) 
a k-point grid of at least 18x18x18 is needed. Taking the av- 
erage over all the £ d values up to a particular NxNxN (dot- 
dashed line) is unhelpful, but taking a weighted average 

helps dramatically, as it effectively increases the k-point den- 
sity: the points in a 4x4x4 grid are not contained in a 6x6x6 
grid, for example. There are two obvious choices for vv/y: the 
best is wn = A/ 3 , the number of k-points in the full Brillouin 
zone, but setting wn equal to the number in the irreducible 
wedge is not bad. (Solid and dashed lines respectively.) All 
subsequent results will be weighted averages using wn = A/ 3 . 
(Incidentally, for the unrelaxed neutral vacancies we find 5 er- 
rors of 0.36 eV and 0.06 eV respectively for Vp and Vi n in 
the 512 atom cell when comparing F point only calculations 
to converged values, so the F point is never sufficient.) 
In reference 11 we also showed that the relaxation energy 

£ Re iax (AO = £ d C:Rx (AO - £ d C:Id (AO (4) 

(where £ d :Rx (A/) and E$ :ld (N) are E$(N) with atoms at re- 
laxed and ideal positions respectively) converges faster with 
k-point grid than £ d :Rx (A/) and £ d :ld (A/) themselves. The re- 
laxed formation energies E^' Rx are then approximated by 

E^ « £p - e Relax (A/) =~E^ + £ d C:Rx (AO - £ d C:Id (N) (5) 

The relaxation energies used will be weighted averages us- 
ing 6x6x6 and 8x8x8 Monkhorst-Pack k-point grids in the 8 
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FIG. 2: Scaling the formation energy for the relaxed and unrelaxed 
structures of the neutral antisite defects Pi n and Inp with inverse su- 
percell size. The x-axis scale is in units of 1/(5. 827 A), the inverse 
of the 8 atom supercell size. Hence x =1.00, 0.50, 0.33 and 0.25 
correspond to the 8, 64, 216 and 512 atom cells respectively. 



atom cell, 2x2x2 and (if the convergence is uncertain) 4x4x4 
grids in the 64 atom cell and 2x2x2 in the 216 and 512 
atom supercells. For the latter two we usually restrict the k- 
point grid to the irreducible Brillouin zone of the undisturbed 
bulk lattice. In other words, we use just the special k-point 
(0.25,0.25,0.25). This amounts to assuming that the distortion 
in the bandstructure due to the presence of the defect is either 
localized (thus important only very near F) or symmetric. It 
introduces a small error whose significance disappears in the 
large supercell limit. 



III. SCALING RESULTS. 

In Figs ElEl and E| we show the formation energies for the 
antisites, phosphorus interstitials and indium interstitials re- 
spectively, both relaxed (minimum energy configurations) and 
unrelaxed (atoms at ideal bulk lattice sites) plotted against in- 
verse supercell size. The solid lines are fits of the four points 
to equation^ The y-axis intersect of each of these fits is the 
scaled formation energy E™ corresponding to the formation 
energy of a single isolated defect in an otherwise perfect, in- 
finite lattice. The inclusion of formation energies from the 
8 atom supercell could be questioned, since in itself it is so 
far from being converged. However, the results shown in the 
figures clearly support our expectation that the correct form 
for the scaling has (at least) three parameters. We therefore 
need results from at least four supercells of the same sym- 
metry. It would have been preferable to use the 1000 atom 
simple cubic cell, but sufficiently converged ab initio calcula- 
tions for InP defects in this cell are not possible with current 
facilities. In addition, the results themselves justify the use of 
the 8 atom cell: the scaling mostly works very well, produc- 
ing small error bars and with E® os lying on or near the curves. 
This, incidentally, also tells us that the k-point convergence in 
the individual cells was sufficient. The cases where the seal- 
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ing does not work so well turn out to be due to other problems, 
see sectionlVl 

To get an idea of the accuracy of the fitting and of the de- 
rived E°f values (and the individual E^) four more fits (dotted 
lines) are added in each case. For each of these, one of the 
four data points has been omitted. The spread in the resulting 
y-axis intersects gives an error bar for the scaled formation en- 
ergy E°£. (It should be emphasized that the dotted lines - three 
parameters to three points - are not intended to be meaningful 
in themselves but are merely indications of the scale of the 
error in the real solid line fits.) The scaled formation energies 
are listed in tableU together with those of the vacancies taken 
fromRef[5]. 

The scaling curves also have a more general and rather 
practical meaning: they are essentially predictions of the for- 
mation energies in all simple cubic supercells from 8 atoms to 
infinity. For example, table HJ shows the predicted formation 
energies in the 1000 and 8000 atom supercells, as they would 
arise from k-point and basis set converged LDA calculations. 
(The 64,000 atom cell would be just as easy, but these predic- 
tions are more likely to be tested and - hopefully - confirmed 
within our lifetimes!) As a more immediate test, table UTI also 
shows the formation energies in the 512 atom cell predicted 
by scaling the results from only the 8, 64 and 216 atom cells: 
following along that dotted line which doesn't pass exactly 
through the 512 atom value for each case in the figures. The 
error in this prediction (as compared to the actual calculated 
values) is also shown. The errors are pleasantly small, espe- 
cially considering only 3 cells have been used, including the 8 
atom one. The errors in the other predictions are expected to 
be much smaller still. 

All the relaxed structures turn out to be symmetric, with 
just breathing mode relaxations. The exception is Inp, which 
shows some moderate Jahn-Teller behaviour (see section lVAl 
) All the interstitials apart from relax outwards, as does 
Inp, whilst Pi n relaxes inwards. This is all as expected, since 
P is smaller than In. Fig[5]shows the scaling of the percentage 
volume change upon relaxation. (The volume shown is that 
of the polyhedron defined by the nearest neighbours.) The fits 
are described in section HVl and the scaled results are in table|I] 
The error bars are derived in the same way as for the formation 
energies, although the dotted lines are omitted here for clarity. 

The most stable neutral native defect in stoiciometric InP 
turns out to be the phosphorus antisite Pi n closely followed by 
the vacancy Vp and then the indium antisite Inp . Then come 
the interstitials, of which indium is the more stable. The least 
stable neutral defect is Vi n . The most stable site for the phos- 
phorus interstitial is the quasi-hexagonal site P;(hex)- For in- 
dium, the two tetrahedral sites are degenerate to within the er- 
ror bars but numerically the phosphorus surrounded site Ini(p) 
is 0.05 eV lower - which is probably correct, since the error 
bars are asymmetric. (See Fig@]) This is a result which is only 
apparent in very large supercells and the degeneracy only ap- 
pears at the 512 atom supercell. In smaller cells In;(i n ) seems 
more stable. What is more, scaling shows that simply taking 
the 512 atom result would give a relaxed formation energy 
about 0.4 eV too high (0.8 eV too high if we stopped at 64 
atoms) which can be large enough to make real differences in 
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FIG. 3: Scaling the formation energy of neutral phosphorus intersti- 



tials. 



Those at the tetrahedral sites Pi(i n ) and P;(p) are not stable in 



the 216 and 512 atom supercells: if the symmetry is not restricted 
then the interstitials relax towards the hexagonal site. Plotted points 
are from relaxations in which the symmetry was restricted to Tj to 
prevent this. 



the predicted physics of InP. As with the formation energy 5 of 
Vin this emphasizes the value of finite size scaling the results 
of supercell calculations, since the largest cells for which we 
can actually do calculations can still be too small to be fully 
converged. 

There are complications with the stable interstitial sites as 
a function of supercell size. Whilst both interstitials are stable 
at all three locations in the two smallest cells they are not so 
in the larger cells. For indium the quasi-hexagonal site lies 
about ~ 1 j eV above the tetrahedral sites and is not stable 
in the larger cells. An indium atom placed here migrates to 
a tetrahedral site upon relaxation, indicating that the quasi- 
hexagonal site is probably not even metastable for indium in 
the real material. For the phosphorus interstitial, on the other 
hand, it is the tetrahedral sites which are unstable. In this case, 
however, we can still obtain (hypothetical) relaxed formation 
energies at the tetrahedral sites by only allowing Td breath- 
ing mode relaxations. The tetrahedral sites are 1.2 eV higher 
than the quasi-hexagonal site. The reason for this difference 
in stabilities may be (partially) stearic: the unrelaxed nearest 
neighbour distances are shorter at the hexagonal site, where 
the smaller phosphorus interstitials are stable and longer at 
the tetrahedral site where the larger indium interstitial sits. 



IV. THE CORRECT FORM FOR THE SCALING 
EQUATION. 

A. Formation and relaxation energies. 

So far it has been assumed that the correct functional form 
for scaling is that of equation ^ This need not be the case, 
however. Equation ^ is based upon approximations and pre- 
dictions for the form of the errors arising from electrostatic 
charge multipole interactions between the defect and it's im- 
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FIG. 4: Scaling the formation energy of neutral indium interstitials. 
Indium is not stable at the hexagonal site Inj(HEX) m tne 216 and 
512 atom supercells and cannot be forced to stay put by restricting 
relaxation symmetry. 
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FIG. 5: Scaling of the volume change with supercell size: percentage 
volume change upon relaxation is plotted against inverse supercell 
size. Pi„ and Pj(p) relax inwards, the others outwards. 



TABLE I: Scaled relaxed and unrelaxed (ideal lattice sites) formation 
energies E^, plus the scaled percentage volume change (%<5V) upon 
relaxation for neutral native defects in InP, listed in order of (relaxed) 
stability in stoiciometric material. (The volume is that of the polyhe- 
dron defined by the nearest neighbours.) Note that the error bars are 
not actually symmetric: see Figs|2|to|4] 



Defect 


Ideal (eV) 


Relaxed (eV) 


%<5V 


Pm 


2.49±0.02 


2.28±0.03 


-19±7 


v P 


3.00±0.10 c 


2.35±0.15 c 


-43±4 C 


lnp 


3.37±0.09 


2.69±0.06 


17±3 


Irii(P) 


4.15±0.30 


2.95±0.20 


19±3 


Pi(HEX) 


4.90±0.14 


3.69±0.08 


12±38 


Pi(P) 


5.21±0.16 


4.88±0.25 a 


-5±7 a 


Pi(In) 


5.49±0.06 


4.96±0.02 a 


9±l a 


in i(In) 


4.75±0.35 


3.00±0.08 


45±19 


I n i(HEX) 


6.95±0.01 


~3.5 b 


~27 b 


Vl„ 


4.95±0.10 c 


4.20±0.05 c 


-35±3 C 



a Unstable in some cells, value results from symmetry restricted 
relaxations, see text. 

b Unstable in some cells, value is rough extension with no error bar 
available, see text. 

c Data taken from Ref. [5] Scaling re-done using n = 2. 



ages in the PBCs and these could be incorrect. In addition, we 
include here relaxations within finite sized supercells, so we 
have both elastic errors and, potentially, cross terms between 
the elastic and electrostatic errors, so other possible scalings 
should be considered. There is clearly a linear term present, 
plus at least one higher order term, so we consider scaling of 
the form: 

£ d c (L) = £ d " + aiZT 1 + a n L- n (6) 

with n = 2, 3 and 4. The fit quality is assessed using the "X 2 " 
test. This is not easy to do reliably with only four points, so we 
average over different defects. Simply summing or averaging 



TABLE II: Predictions for the unrelaxed and relaxed formation ener- 
gies in the (to date) uncalculated 8000 and 1000 atom simple cubic 
supercells. Also shown are predictions for the 512 atom supercell 
formation energies obtained by scaling the values in the 8, 64 and 
216 atom cells, together with the difference (Error) between the pred- 
icated and calculated values. (Energies in eV.) 



Idea Structures Relaxed Structures 



Defect 


8000 


1000 


512 


Error 


8000 


1000 


512 Error 


Pm 


2.50 


2.51 


2.51 


-0.01 


2.31 


2.33 


2.35 0.00 


v P 


2.99 


2.98 


3.00 0.04 


2.40 


2.46 


2.46 -0.03 


lnp 


3.29 


3.21 


3.18 


0.02 


2.68 


2.66 


2.66 0.01 


Ini(P) 


4.30 


4.44 


4.55 


0.07 


3.13 


3.31 


3.43 0.05 


Pi(HEX) 


4.90 


4.89 


4.91 


0.03 


3.68 


3.67 


3.68 0.02 


Pi(P) 


5.13 


5.05 


5.05 


0.06 


4.81 


4.74 


4.75 0.06 


Pi(In) 


5.49 


5.48 


5.45 


0.02 


4.99 


5.00 


5.00 0.00 


I n i(In) 


4.87 


4.97 


5.08 


0.08 


3.15 


3.31 


3.40 0.02 


I n i(HEX) 


6.98 


7.01 


7.02 


0.00 








Vm 


4.88 


4.80 


4.77 


0.02 


4.19 


4.18 


4.17 -0.02 



the x 2 values would bias the conclusion towards the worst 
data sets. Instead, for each data set we find x 2 f° r eacn value 
of n, select the n giving the best fit Gtu J an d tnen calculate 
a quality factor Q n =X 2 ^X^, est f° r eacn n - We then compare the 
averaged Q n . 

We first examine the scaling of the elastic errors. To do this 
we have performed a series of relaxations in the 216 atom cell 
for each of six defects. In these relaxations, the number of 
shells of atoms permitted to relax around the defect is varied 
from 1— >4 for interstitials or 1— >5 for antisites and vacancies. 
(4 and 5 are the maximum numbers of complete shells which 
fit inside the supercell.) Since the cell size is kept constant 
the electrostatic interactions will be (almost) constant so any 
variation in E& is due to elastic effects. In the left panel of 
Fig |6] we show the scaling of the relaxation energy with the 
inverse radius of the outermost relaxing shell. (The units have 
been scaled to match those in previous figures.) The scaling is 
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TABLE III: Scaled percentage volume changes upon relaxation, 
from scaling with number of shells relaxed in the 216 atom cell and 
from scaling with supercell size. 



Defect Shells Superells 



Pin 


-15±1 


-19±7 


Illp 


20±2 


17±3 


Ini(p) 


17±1 


19±3 


I n i(In) 


51±40 45±19 


Pi(HEX) 


23±2 


12±38 


V P 


-39±3 


-43±4 



almost purely linear (solid lines in the figure) even - within the 
bounds of error - for the Jahn-Teller active defects like Inp and 
Vp. Adding a higher order term, such as L~ 3 (dashed lines) 
clearly only improves the fit very slightly: % 2 is reduced by 
about 30% on average. 

Having established that the scaling of the relaxation energy, 
and hence of the elastic errors in £^°, is linear we turn to the 
scaling of the formation energies with supercell size. Unfor- 
tunately there is too much scatter in several of the curves so 
that whilst n = 3 is best in some cases, in others it comes an 
often close second to n = 2 or n = 4. However, since the elas- 
tic errors are now known to be linear in relaxation radius and 
hence in supercell size, we can assume that it will simply add 
to the linear term in equation |6] so that we can calculate Q„ 
using all 19 data sets (Figs [2] to [4] and Ref. [5]). The result 
obtained is that the scaling does indeed fit best with n = 3, 
with both Q2 and Q4 being 4 times larger than Q3. In fact, 
if we calculate the {Q„} using only the unrelaxed formation 
energies we still find n = 3 fits best and if we use the relaxed 
formation energies minus that of Pr n we find the same result. 
Unfortunately if we use all of the relaxed energies including 
Pin we get a different result: n — 2 provides a better fit. This 
suggests the (faint) possibility of a cubic scaling for unrelaxed 
formation energies but quadratic for relaxed energies, which 
shifts the predicted by 0.01^0.2 eV depending on the de- 
fect. A quadratic scaling would seem odd, since it comes from 
neither the elastic nor the electrostatic errors but it could pos- 
sibly arise from the cross terms between them. An alternative 
and perhaps more likely explanation is that the uncertainty re- 
garding n is a result of the spurious dispersion of the defect 
levels mentioned above. This adds an exponential term to the 
formation energies in the smallest supercells, blurring the pic- 
ture slightly. It seems most likely that the scaling should be 
cubic even for the relaxed formation energies: even now it is 
only one defect which seems to particularly disagree. This 
should be confirmed once reliable calculations involving the 
1000 atom supercell become feasible. 



B. Relaxed volumes. 

Fig [5] shows the scaling of the percentage volume change 
(going from ideal to relaxed structures) with supercell size, 
whilst the right panel of Fig [6] shows the volume changes 




1/shell radius (1/8 atom cell size.) 



FIG. 6: The relaxation energy (£R e i ax ) plus the % volume change 
as a function of the number of shells around a defect which are al- 
lowed to relax. For Inpthe % Jahn-Teller distortion in the ADX and 
DAD structures (see main text) is also shown. (Relaxation energy 
and volume change are only shown for the DAD structure: those for 
the ADX structure are very similar.) The 216 atom supercell is used, 
giving 1— *5 shells for vacancies and antisites, 1— >4 for interstitials. 
The x-axis is the inverse of the radius of the atom shell, in units of 
one over the 8 atom cell size. Volume changes are outwards for all 
except P[ n . Fits are to equation [6] with (dashed) and without (solid) 
the L~" term, (n = 2 for £R e iaxi » = 3 for the volumes.) Results for 
Vp are taken from Re». 



scaled with the number of shells relaxed in the 216 atom su- 
percell. For the antisites and In;(i n ) a linear fit is again rather 
good, but for all of the other data sets a higher order term 
is clearly present. We have again tried adding terms in L~ n , 
n = 2, 3 and 4. For both supercell and shell scaling a L~ 2 
is best. For supercell scaling Qt, is 42 times larger than Q2 
and Q4 is 124 times larger. The curves plotted in Fig [5] are 
thus quadratics, with the y-axis intercepts giving the volume 
change expected for a lone defect in an infinite supercell. 
These scaled volumes are shown in table U For shell scal- 
ing (Fig[6} Q3 is 1.4 times larger than Q2 and Q4 is 2.7 times 
larger. The quadratic fittings are shown as dashed lines in the 
right hand panels of[6] The y-axis intercepts this time give the 
volume change expected if an infinite number of shells were 
relaxed. They are shown in table|lll] where the equivalent val- 
ues for supercell size scaling are added for comparison. The 
two sets of values agree very well, indicating that at least the 
breathing modes in the infinite limit are unaffected by charge 
multipole interactions. 



V. WHY SOME STATES SCALE BETTER THAN OTHERS. 
A. Jahn-Teller active defects. 

It is very clear from table H] that some formation energies 
scale better than others. When considering the scaling for the 
neutral vacancies 5 we noted that the scaling becomes more 
difficult to do reliably for strongly Jahn-Teller active defects 
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such as Vp. For a Jahn-Teller active defect there is a par- 
tially filled degenerate state at the Fermi level, which the Jahn- 
Teller theorem 6 says will be lifted by symmetry reducing re- 
laxations (if no other effect achieves this first). This leads to 
poor scaling since the symmetry of the most stable relaxed 
structure can vary with supercell size, so that data points from 
some cells scale differently to those from others. In order to 
get good error bars for scaled formation energies each pos- 
sible reduced symmetry structure must be scaled separately. 
Amongst the current defects there is a further example of a 
Jahn-Teller active defect, In P . The distortions here are much 
weaker, so the error bar is ±0.06 eV even when symmetry dif- 
ferences are ignored, which is still reasonable. Never-the-less, 
we have done a search for the various stable and metastable 
structures in the four supercells. Their various formation en- 
ergies are shown scaled in FigQ together with a) the scaling 
of the lowest lying formation energy irrespective of symme- 
try (labelled Min) and b) the formation energy when only T c / 
symmetry breathing mode relaxations are permitted (labelled 
Td). 

In the 8 atom cell the lowest lying structure has the full T ( / 
symmetry of the unrelaxed antisite. Indeed, this is the only 
stable structure we find in this cell. In the larger cells relax- 
ation breaks the T c i symmetry to give (primarily) C3,,, T>2h or 
C2 V point groups at the defect site. C3 V symmetry is reached 
by the antisite atom moving either towards the midpoint of 
three of its nearest neighbours (away from the fourth) in a DX 
like structure (DX in the figure), or towards one neighbour and 
away from the other three - an anti-DX (ADX) structure. Q2 V 
arises when the antisite moves towards one pair of neighbours 
and away from the other pair, i.e. towards a bond centre site 
(BCR). The degeneracy can also be lifted when the antisite 
stays still with the four antisite-neighbour distances equal, but 
the angles between the bonds changes. D2/, structures occur 
if the neighbours rotate to form two opposing pairs of either 
shorter or longer neighbour-neighbour distances (DDM - dou- 
ble dimer or DAD - double antidimer, respectively). In the 64 
and 512 atom cells the lowest lying structure is a DAD struc- 
ture with two neighbour-neighbour distances respectively 5% 
and 8% shorter than the others. The most stable structure in 
the 216 cell is a 7% DAD-like structures but with an addi- 
tional 4% DX-like distortion, although a 4% BCR structure 
and a 7% pure DAD structure (with no DX component) both 
come a close second. (The distortion quoted for the BCR, DX 
and ADX structures is the % variation in antisite-neighbour 
distances.) In the 512 atom cell the potential energy surface 
for small (up to ~2%) DX distortions from the DAD struc- 
ture is also very flat. Overall, these results suggest that a lone 
Inp in an infinite supercell would have a DAD structure with 
a formation energy lying about 0.4 eV below the Tj structure 
found when only breathing mode relaxations are allowed, and 
0.1 eV below the formation energy found by scaling the min- 
imum formation energy irrespective of Jahn-Teller structure. 

The changes in relative stability of the different structures 
are due to one or a combination of two things: a) stabiliz- 
ing/destabilizing dipolar, quadrupolar or higher interactions, 
which can in certain cases lift the symmetry without distor- 
tion, (in the 8 atom cell for example,) or favour certain Jahn- 
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FIG. 7: Scaling the formation energies for various Jahn-Teller struc- 
tures of Inp. In all panels the global minimum is shown (Min). In 
addition we show the formation energy when only breathing mode 
relaxations are allowed (Td) and non-symmetry restricted relaxed 
formation energies for the BCR, DX, ADX and DDM structures, 
plus the DAD structure with (DAD-DX) and without (DAD-EQ) an 
additional DX like distortion. (See main text for descriptions.) Fits 
to equation \\\ are solid when the structure is stable in four cells or 
dashed when it is only stable in three. Dotted lines are linear fits 
when the structure is stable in only two cells. 



Teller structures over others. These effects become weaker as 
the cells grow, b) The lack of shells of atoms in the smaller 
cells to absorb the elastic strain, which favours more symmet- 
ric structures. In the right hand panel of Fig [6] the variation 
in the degree of Jahn-Teller distortion for the ADX and DAD 
structures was plotted versus the number of shells permitted 
to relax within the 216 atom cell. For the ADX structure there 
is virtually no variation at all and the same was previously 5 
found for Vp. For the DAD structure, on the other hand, a 
rather strong variation is found. This suggests that elastic ef- 
fects are involved for some local distortion symmetries but not 
for others, thus further complicating the possible variations in 
lowest symmetry structure with supercell size. At least for the 
ADX structure the charge multipolar interactions act essen- 
tially in competition with the normal Jahn-Teller mechanism, 
making it uncertain if the correct structure has been found at 
all for smaller supercells. It is sometimes pointed out (Ref 
[14] and elsewhere) that one way around this problem would 
be to use k-point integration at the F point only, since this 
restores the degeneracy of the degenerate levels (prior to re- 
laxation). However, since for InP (and doubtless many other 
materials) the F point does not give sufficiently converged 
formation energies even in the 512 atom supercell this will 
simply result in unconverged results: errors arising from the 
use of just the F point can be tenths of eV, often the same 
size or larger than the splittings between various stable and 
metastable Jahn-Teller distorted structures. Here we study the 
relative stability of all the possible Jahn-Teller distorted struc- 
tures with converged k-point grids and then scale them to pre- 
dict which structure will be most stable in the infinite super- 
cell limit, where these spurious degeneracy lifting interactions 
become zero. 
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FIG. 8: Defect levels in or near the band gap, shown at the T point for neutral native defects in the 64 and 512 atom supercells. Dotted 
horizontal lines mark the bulk LDA valence and conduction band edges. (0.6771 eV apart.) Dashed lines indicate the shifted bulk and valence 
band edges in the defect cells, solid lines indicate localized defect levels. Black, grey and white symbols are (respectively) filled, half filled 
and empty states. Circles indicate localized states, diamonds delocalized (band) states. Filled valence band and empty conduction band states 
are not marked with symbols. 



B. Defect level dispersion and defect states outside the 
bandgap. 

The scaling is also rather bad for the interstitials, even 
though no Jahn-Teller behaviour is anticipated or detected 
what-so-ever. One might have expected that the reason for 
this was the defect level dispersion which was also identified 
in section IIVI as a possible source of the slight uncertainty 
in the correct form for the scaling equation. The dispersion 
can lead to errors in the formation energies of individual de- 
fects and since it should decrease exponentially with cell size 
it could affect the error bars on the scaled infinite supercell 
formation energies. In a fully occupied defect level the dis- 
persion should have no direct effect as long as the defect level 
lies within the gap at all k-points, since the average energy of 
the band should equal the energy of the level in the absence of 
dispersion. On the other hand, the dispersion can lead to hy- 
bridization of the defect level with conduction band states of 
the same symmetry, thus artificially lowering the mean value 
of the defect level and hence of the defect formation energies 
in the smaller cells. (The same holds for empty defect levels 
hybridizing with valence band states, allowing the latter to be 
lowered in energy.) For partially filled defect levels the effect 
upon the formation energy is more direct, since only the lower 
parts of the defect level (band) will be filled, again leading to 
too low a value for the formation energy. However, the con- 
nection between the amount of defect level dispersion and the 
size of the errors in the formation energy in a particular su- 
percell is not simple. Indeed, checks on the bandwidth of the 
defect levels find no correlation at all between them and the 
scaling error bars: for example the unrelaxed structures for 
both Pi n and In;(i n ) have a defect level in the lower part of the 
gap with a rather large dispersion (^0.6 eV in the 64 atom su- 
percell) and yet the scaling error bar for Pi n is 0.02 eV whilst 



that for In;(i n ) is 0.35 eV. 

The main difference is actually that Inj(r n ) also has a par- 
tially filled defect level inside the conduction band, resulting 
in an electron occupying a delocalized state at the conduc- 
tion band edge. This occurs for all of the defects which have 
poor scaling error bars. This is seen clearly in Fig|8]where we 
present the level diagrams for all the neutral native defects, 
shown at the T point in the 64 and 512 atom cells. While 
the antisites have no delocalized states, the interstitials always 
have at least one electron occupying a state at the conduc- 
tion band edge. For phosphorus at the tetragonal sites not one 
but three electrons lie above the conduction band edge, which 
may be why they are not stable at all in the larger supercells. 
The level diagrams for the vacancies seem at first glance to 
contradict the rule: Vp has one delocalized electron in the 
smaller cells and Vi n has three delocalized holes in all four 
cells, despite the fact that the scaling is rather good for these 
defects. This is itself a finite size effect, however. For Vp 
a triply degenerate defect level lies just above the conduction 
band edge. This level Jahn-Teller splits, with one state moving 
downwards. In the smaller cells it does not make it to the band 
gap so the metastable distorted structures found must arise 
from hybridization between the localized Jahn-Teller split lev- 
els and levels at the conduction band edge. For the largest 
cells, however, the lower Jahn-Teller split level drops into the 
gap, so that the charge state becomes stable. In the case of Vr n 
we again see a triply degenerate localized defect level, this 
time inside the valence band leaving three holes at the valence 
band edge, and even in the 512 atom cell this is all we see. 
However, as the supercell size grows the defect level moves 
rapidly towards the valence band edge. Clearly, for an iso- 
lated defect in a real material we would expect this level to 
lie inside the band gap, even though the 512 atom cell is not 
sufficiently large to show it. This again underlines the need 
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a) 



b) 



d) 



FIG. 9: Schematic diagrams for the electrostatic origin of linear for- 
mation energy error scaling, a) A positively charged defect, with a 
positive core, localized electrons in a defect level and compensating 
jellium. b) Subsystem leading to linear errors: weaker effective core 
charge plus jellium. c) The effective situation for In;(p): core charge, 
two localized defect level electrons and one delocalized electron at 
the conduction band edge, d) The effective situation for the unre- 
laxed Inp: core charge, four localized defect level electrons and two 
(partially) delocalized defect level electrons. 



for using comparison and scaling rather than just large super- 
cells. (Similar effects are visible for Pj (In ) and Pi (P ), although 
they remain unstable.) 



VI. POSSIBLE CHARGE STATES ACROSS THE BAND 
GAP. 

Although we only consider the neutral defects in this pa- 
per, the level diagrams in Fig [8] also give a rough indication 
of which charge states are most likely in different parts of 
the band gap. For the indium interstitials the +1 charge state 
seems the most stable in the upper part of the band gap, but 
two transition levels will lie in the mid-gap, giving a +3 charge 
state for p-type material. For Pi(hex) we expect a +1 charge 
state across most of the gap. With increasing supercell size 
a second defect level approaches the band edge from below 
so there may be one or two transition levels near the bottom 
of the gap. This is a results which is only apparent when one 
compares or scales the results from different supercell sizes: 
even in the 216 or 512 atom cells it is not apparent. Pi n has 
a single, filled level in the upper half of the gap, suggesting 
for n-type material and +1 or +2 otherwise. Inp has a three- 
fold degenerate level containing 4 electrons in the middle of 
the gap. A large number of transition levels across the entire 
gap are thus expected, ranging from -1 or -2 in n-type material 
to potentially even +4 for strongly p-type material. Vp should 
have two transition levels in the upper half of the band gap, 
+ 1 and -1 being the most stable charge states in strongly p and 
n-type material respectively. For Vr n we expect the -3 charge 
state to be stable from the mid-gap upwards. However, the 
movement of the threefold degenerate defect level up into the 
band gap (section lV Bi leads us to expect six transition levels 
all lying near (above or below) the valence band edge. The 
most stable charge state at the valence band edge itself could 
be anything from to +3. 



VII. THE ORIGIN OF THE LINEAR SCALING TERM 
FOR NEUTRAL DEFECTS. 

We expect linear scaling terms to arise from two sources: 
from incomplete relaxation (elastic errors) and from electro- 
static errors. For electrostatic errors linear scaling comes 2 
from the monopole term in the multipole expansion of 
the electrostatic interactions: roughly speaking, it is the 
Madelung energy of the localized charge on the defect inter- 
acting with the compensating jellium background in the image 
cells. This is shown schematically in Fig [9] Part a) shows the 
typical situation for a positively charged defect: An infinite 
array of very tightly localized positive core charges (one in 
each periodic image of the supercell) are only partially com- 
pensated by the localized electrons in the (less tightly) local- 
ized defect levels. The remaining charge is compensated by 
jellium. Makov and Payne 2 extracted the part of the system 
shown in Fig|9]b) consisting of an infinite array of positive 
delta functions interacting with jellium and showed that is 
gives rise to a formation energy error which is linear in the 
supercell size, with a strength proportional to the square of 
the charge on the defect. This error should therefore be ab- 
sent for unrelaxed neutral defects. The fact that we see a clear 
linear contribution to the unrelaxed formation energy of, for 
example, the neutral Inp and tetrahedral In interstitials is thus 
somewhat unexpected. 

The explanation lies in the localization/delocalization of the 
states which become occupied and unoccupied when the cell 
contains a neutral defect. In the cases of the neutral intersti- 
tials and vacancies one or more defect level(s) lie outside the 
bandgap for all or some supercell sizes. As a result there are 
either electrons in delocalized band states at the LDA con- 
duction band edge or holes at the valence band edge. This 
is shown for the example of In^p) in Fig ^| To quantify the 
localization of a particular Kohn-Sham eigenstate we start by 
projecting it onto Wigner-Seitz cells 15 around each atom. We 
then select some radius rj around the defect and find the av- 
erage projection per atom p> over the part of the cell where 
r > rj, and similarly the average p< for the r < rj region. The 
localization in then given by the ratio 



a, 



(l) 



P> 
P< 



(7) 



For an isolated defect in an otherwise perfect infinite lattice, 
CCp — > oo as rj — > °° if the state is localized, since p> -^0. On 



the other hand, a n 



1 if the state is delocalized as the two 



averages p> and p< then tend to the same value. 

Here we have finite L sized supercells, so we choose rj 
such that p> is calculated over the outermost complete shell 
of atoms around the defect, plus the atoms in the cell corners. 
p< is then calculated over the remaining 2, 4 and 6 complete 
shells in the 64, 216 and 512 atom cells respectively. (The 
8 atom cell is too small for this analysis.) Supercell scaling 
of Otp(rj(L)j then has the same behaviour as for the isolated 
defect. Thus, in Fig a p is plotted against one over the 
cell size for various states of interest. The left panel shows 
the localization scaling of the states at the LDA valence band 
maximum (VBM) and conduction band minimum (CBM) for 
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1/L (1/8 atom cell size) 



FIG. 10: Localization of the Kohn-Sham levels corresponding to 
(left panel) the valence and conduction band edges at the T point in 
the cell containing In^pj; (middle panel) the defect levels in the band 
gap for Inj(p), Pi n and Inp, calculated at the F point; (right panel) 
the three Inp defect levels averaged over the whole Brillouin zone. 
See text for details. Lines are simple quadratic fits for guidance: the 
correct form for this scaling has not been investigated. 



cells containing an unrelaxed In;(p), whilst the middle panel 
(with a very different vertical scale) shows it for the mid-gap 
defect levels of both Inj (P ) and the antisites. The well behaved 
Pi n is included for comparison. These localizations, like the 
level diagrams in Fig [8] have been calculated at the F point, 
but using fully k-point converged charge densities. For In;(p) 
the band states are delocalized and the defect level is local- 
ized. The VBM and defect levels are filled, whilst the delo- 
calized CBM level is half filled. Hence adding the defect to 
the cell has added 2 localized electrons and one delocalized 
electron in the vicinity of the band gap. Hence, we have an 
electrostatic situation like that shown schematically in Fig [9] 
c). To a first approximation we can replace the charge density 
of the delocalized electron by its average value, thus recov- 
ering the situation of Fig|9]a). We thus predict a linear term 
in the formation energy scaling. However, the distribution of 
the delocalized charge is in reality far from uniform on the 
scale of the atomic spacing, so it is not clear what the prefac- 
tor and effective charge should be. Hence, without performing 
the detailed mathematical derivation required (which lies be- 
yond the scope of the current paper) we cannot predict what 
gradient the linear term should have. 

The case of Inp is different. There are no electrons in de- 
localized conduction band states or holes in the valence band. 
However, the localization of the defect level at the F point is 
rather weak, certainly compared to Pi n and In; ( p) . The defect 
level of Inp is three-fold degenerate and partially filled. Away 
from the F point this degeneracy is split by the interaction of 
the defect with its images, leading to different dispersions of 
the defect levels in different directions in k-space, with one 
level having lower energy (away from F). At most k-points 
this state is thus completely filled, whilst the other two are 
partially empty. Hence, averaging the Wigner-Seitz projection 
over the whole Brillouin zone (using 8x8x8, 4x4x4 or 2x2x2 
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FIG. 1 1 : Scaling of errors from other sources. Left panel: the change 
in formation energy when the In 4d electrons are treated as core 
rather than valence - i.e. the formation energy using an US-PP with 
In 4d treated as valence minus the formation energy found earlier 
with the In 4d as core. Centre panel: the change in formation energy 
when PAW is used rather than US-PP (both with In 4d are treated as 
valence.) Right panel: the change in formation energy when GGA is 
used instead of LDA. Solid line fits are to Eq|8| dashed lines are to 



point Monkhorst-Pack grids in the 64, 216 and 512 atom cells 
respectively) the completely filled level becomes fully local- 
ized. (Right panel of the figure.) The other two levels are more 
occupied at the origin in k-space than elsewhere and hence in 
real space are only partially localized. Clearly the Jahn-Teller 
structural relaxations are required in order to properly local- 
ize these states. In the unrelaxed structure we effectively have 
two localized electrons and two partially delocalized ones in 
the middle of the band gap, leading to the spatial charge distri- 
bution shown schematically in Fig[9jl). Electrostatically, this 
again behaves to first order like a jellium background charge 
surrounding a positively charged, localized defect, resulting 
in the linear scaling observed. As with the interstitials, pre- 
dicting the specific gradient expected must be left to future 
work. Here, we simply note that this mechanism should also 
affect the unrelaxed formation energies of other defects with 
partially filled degenerate levels in the band gap. It seems rea- 
sonable to presume that it is involved with the neutral Vi n and 
perhaps Vp, although they also have linear terms arising from 
partially filled band states. 



VIII. OTHER SOURCES OF ERROR. 

Finite size errors are not the only types of errors present in 
these calculations: errors also arise from the truncation of the 
basis set and the k-point integration as well as from the use of 
both pseudopotentials and LDA. Furthermore, memory size 
limitations forced us to use pseudopotentials in which the in- 
dium 4d shell is treated as core rather than valence. Although 
the central aim of this paper is to study the treatment of fi- 
nite size errors it is still informative to estimate the size of 
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these other sorts of errors for comparison. In principle, errors 
arising from the pseudopotentials and the LDA should be in- 
dependent of supercell size, though some short range depen- 
dence may still be anticipated since changes may affect the 
amount of defect band dispersion. This contribution should 
disappear exponentially with supercell size. 

In fig^2 we show the change in the unrelaxed formation en- 
ergies when (left panel) the In 4d electrons are treated as core 
rather than valence, (centre panel) the US-PP are replaced 
by the Projector Augmented Wave method 16 (PAW) or (right 
panel) the LDA exchange correlation functional is replaced 
by the Generalized Gradient Approximation (GGA) of Perdue 
and Wang 17 . These formation energy differences are shown as 
a function of supercell size for the seven stable point defects. 
(Vin, In P , Ini(i n ), Ini(P), Vp, Pi n and Pi (H EX)-) For the GGA and 
for LDA with the In 4d as valence a planewave cutoff of 200 
eV gives the same level of accuracy as it does with the In 4d 
as core. For PAW on the other hand a planewave cutoff of 
300 eV is required for the same accuracy level. The relaxed 
lattice constants and bandgaps also change. For LDA with In 
4d as valence we find 5.833 A and 0.581 eV respectively, for 
LDA with PAW we find 5.830 A and 0.597 eV whilst for GGA 
(with US-PP and In 4d as core) we find 5.956 A and 0.473 eV. 

The errors coming from the pseudopotentials (left and cen- 
tre panels of Fie ll 1> are, as anticipated largely independent of 
the supercell size. The slight dependence is well described by 
a two parameter fit to the general exponential 

£ d c (L) = E" + a e (exp (IT 1 ) - 1 ) (8) 

(Solid lines in the figure.) Changing the In 4d electrons from 
core to valence has, not surprisingly, only a fairly small effect 
(well below 0. 1 e V) upon the phosphorus related defects, but 
a much more significant effect upon the indium related ones, 
particularly In;(i n ) and Inp where the effect is on the same or- 
der of magnitude as the finite size errors. PAW produces more 
accurate results than any pseudopotential method since it re- 
constructs the exact valence wavefunction with all nodes in 
the core region. The replacement of US-PP with PAW pro- 
duces only small (O(0. 1 eV) or less) changes in the formation 
energies with virtually no size dependence: the largest differ- 
ence between the correction in the 8 atom cell and that in the 
64 atom cell is only 0.01 eV, so the US-PP to PAW changes 
have not been calculated in the 216 atom cell. This demon- 
strates that the widely used US-PP 9 are perfectly reasonable 
for this type of calculation. It should also be noted that in most 
cases the small correction when PAW is introduced partially 
cancel that made when moving the In 4d electrons from core 
to valence, at least for the examples studied here. 

The errors arising from the exchange-correlation functional 
have two main forms. Firstly, the band gap is (usually) 
strongly reduced compared to experiment. (For InP GGA 
does worse than LDA once lattice parameter optimization has 
been included.) This reduction leads to ambiguities in the def- 
inition of the formation energy for charged defects, in turn 
leading to large uncertainties in predictions. For some semi- 
conductors the bandgap can even be reduced to zero mak- 
ing the material appear metallic and dramatically altering the 
properties of many defects. However, neither of these ef- 



fects occurs here, since we consider neutral defects and ob- 
tain a non-zero band gap. Never-the-less, a second exchange- 
correlation related error is present: LDA overbinds all bonds, 
moving some defect formation energies up and others down. 
Hence to assess the errors involved we here compare with 
GGA, which is known to have the opposite effect: underbind- 
ing. An exact solution to the DFT formation energies would 
lie somewhere in between - and probably closer to the LDA 
results since LDA gives a better lattice constant for bulk InP. 
We note also that we have not allowed spin polarization in 
our calculations, using LDA instead of the Local Spin Den- 
sity Approximation (LSDA). However, for most of our defects 
the use of LSDA would simply cause further finite size errors 
since it would introduce spurious magnetic interactions be- 
tween the defect and it's PBC images. The exceptions occur 
for the Jahn-Teller active defects (Inp, Vp and Vi n ) where the 
degeneracies can in certain cases be lifted by Hund's rule cou- 
plings. However, this is only important for materials 18 with 
more tightly localized bonds or dangling bonds than we have 
here, so we may safely omit it. 

The results from replacing LDA with GGA are shown in the 
right panel of Fig^2 an d arc somewhat surprising: there is a 
clear finite size error associated with the choice of exchange- 
correlation functional for certain defects, namely In^p), Ini(i n ), 
Vp and Pi(HEX)- The antisites on the other hand show no sig- 
nificant supercell size dependence. (Vi n shows only a small 
size dependence but is still fitted better by a polynomial than 
by an exponential.) The defects for which the LDA/GGA er- 
ror difference depends on supercell size are those for which 
the defect levels lie outside the band gap. This indicates that 
the spurious delocalized states are treated differently by dif- 
ferent exchange correlation functionals. This is confirmed by 
checking the degree of localization of the electrons/holes at 
the band edges: when GGA is used the values of a p are 20 
to 400% larger than with LDA. These changes are linked to 
the change in band gap and also to the fact that the defect 
bands move closer to the band edges (at least near the F point) 
for these particular defects. Hence we may expect large, su- 
percell size dependent errors of 0(0.5 eV) from the choice 
of exchange correlation functional when partially filled defect 
bands lie outside the bandgap but we expect smaller O(0. 1 eV) 
size independent errors when the defect levels lie within the 
band gap. 

Table IIVI compares the size of the errors arising from dif- 
ferent sources. The basis set errors are estimated from the 
difference in formation energy in the 64 atom supercell when 
the planewave cutoff is raised from 200 eV to 400 eV. The k- 
point errors are different in each supercell: the errors shown 
in the table are the largest occurring for any of the supercells 
used, estimated as the difference between the mean value E^ 
and the value £?(iVmax) obtained using the largest k-point 
grid actually calculated for the defect and supercell in ques- 
tion. All other sources of error - such as defect band disper- 
sion etc - are contained within the finite size errors shown. 
Most of the errors listed are on the 0.1 eV scale or below, 
which in practical calculations is usually acceptable. The fi- 
nite size errors and some of those arising from treating the In 
4d electrons as core are larger, lying on the ^0.5 eV level. 
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TABLE IV: Comparing the size of the errors arising from the vari- 
ous different approximations, i) Finite size errors from the supercell 
approximation (shown for the 64 and 512 atom supercells.) ii) The 
treatment of the In 4d electrons as core (scaled to an infinite cell.) iii) 
US: the use of ultrasoft pseudopotentials, compared to PAW (scaled 
to an infinite cell.) iv) The LDA, compared to GGA, (scaled to an 
infinite cell. For the cases in which a finite size term appears in the 
LDA vs GGA error a more valuable comparison is of the errors in 
individual cells, so the values in the 64 and 512 atom cells are then 
given in brackets.) v) Basis set truncation (in the 64 atom cell) and 
vi) k-point integral truncation. (Shown for the cell in which it is 
worst for the defect in question.) vii) "Final" shows the scaled LDA 
formation energy when the pseudopotential errors are accounted for. 
(All energies in eV.) 



Defect Supercell In4d US LDA Basis K-grid Final 

0.40/0.20 0.11 0.00 0.25(0.24/0.23) 0.003 0.003 4.84 

In P 0.41/0.21 0.30 0.10 0.12 0.003 0.0007 3.57 

In i(In) 0.40/0.25 0.08 0.16 0.48(0.18/0.32) 0.01 0.007 4.67 

In i(P) 0.55/0.33 0.47 0.06 0.32(0.05/0.17) 0.02 0.007 4.56 

Pi„ 0.03/0.03 0.02 0.01 0.07 0.01 0.002 2.46 

V P 0.12/0.03 0.07 0.00 0.40(0.23/0.31) 0.004 0.001 2.93 

Pi(HEX) 0.06/0.02 0.02 0.07 0.13(0.07/0.10) 0.04 0.009 4.85 



We note, however, that for charged defects we anticipate even 
larger finite size errors - up to 1-2 eV in many cases 12 - whilst 
we do not anticipate the errors from the In pseudopotentials 
being significantly different from those here. The presence of 
defect level dispersion effects is confirmed by the existence of 
the exponentially shrinking supercell size dependence in the 
errors related to pseudopotentials and to the LDA for the an- 
tisites. These results also confirm that they are short ranged. 
Their energy scale is rather hard to estimate directly, but the 
size of the exponential components in Fig^2 su gg ests that the 
errors involved are only on the 0.01-0.1 eV scale even for the 
8 atom cell. An alternative estimation comes from the error 
bars on the scalings themselves (see table ^ which are 0.01- 
0.35 eV, only a small part of which comes from the defect 
level dispersion. 

The fact that the errors coming from the pseudopotentials 
have only a small and exponentially decaying cell size de- 
pendence means that it is perfectly reasonable to make the 
approximations we needed to make in order to be able to do 
calculations in sufficiently large supercells to correctly assess 
the finite size errors. The amplitude and sign of these non- 
size dependent errors can be calculated separately in a smaller 
cell - the 64 atom cell for example - and simply added or sub- 
tracted from the final finite size scaled results in order to pro- 
duce much more accurate and reliable defect formation ener- 
gies than those normally published. This has been done for 
the LDA formation energies of the present examples in the 
last column of table IIVI but is equally valid for the errors in 
the relaxed formation energies and those in the GGA (at least 
if no delocalized hole/electron states have appeared at the va- 
lence/conduction band edges.) 



IX. CONCLUSIONS. 

We have studied the finite size errors which occur when 
the supercell approximation is used in the calculation of the 
formation energies and structures of point defects in semicon- 
ductors, using the neutral native defects of InP as an example. 
We have calculated the relaxed and unrelaxed formation en- 
ergies using planewave ab initio DFT in simple cubic super- 
cells containing 8, 64, 216 and 512 atoms - the largest cur- 
rently computable. To examine and correct for these errors 
we have used finite size scaling with inverse supercell size, 
which we consider to be the most reliable and accurate way 
to treat the finite size supercell approximation errors. Unlike 
other methods this does not rely upon any modelling or as- 
sumptions about the errors, other than that they are primar- 
ily long ranged (polynomial rather than exponential) and de- 
crease with increasing supercell size. This method requires 
the results of calculations in at least 3 supercells, and at least 
4 if we are to have an idea of the accuracy of the resulting 
scaled energies. Hence for some difficult cases in which the 
8 atom cell is simply too unreliable it may occasionally be 
necessary to use supercells with up to 1000 atoms. 

Three sources of finite size error have been examined: in the 
case of relaxed formation energies there are elastic errors due 
to the finite volume available for relaxation. We showed that, 
as they should, these scale linearly with inverse supercell size 
(L _1 ) with very little hint of any higher order error term aris- 
ing, even when Jahn-Teller distortions are taken into account. 
The second type of error is the dispersion of the defect levels, 
which has only a relatively small effect upon the formation 
energies, at least when the defect levels are completely filled. 
These effects appear to shrink exponentially with increasing 
supercell size, as anticipated, but appear to slightly increase 
the uncertainties in the final scaled formation energies. The 
third source of error is much more significant and arises for 
both relaxed and unrelaxed formation energies. It is due to 
charged multipole interactions between a defect and its im- 
ages in the PBCs. We have shown that these errors are present 
and not always negligible, even for neutral defects. Linear er- 
rors can arise if fully or partially filled defect levels lie outside 
the band gap in the neutral charge state, leading to delocal- 
ized holes at the valence band edge or electrons at the con- 
duction band edge. Linear errors can also appear in unrelaxed 
formation energies due to the way in which the defect/image 
interactions lift degeneracies for partially filled degenerate de- 
fect states within the band gap. Both of these ways for linear 
scaling errors to arise can apply for other non-neutral defects. 
They indicate that the calculation is in some way unphysical, 
for example that the charge state involved is not actually sta- 
ble. However, in practise it is often necessary to calculate 
formation energies of such unstable charge states in order to 
check which transitions levels do actually lie inside the gap. It 
is thus important to be aware that large finite size errors, such 
as those reported here, can occur even in calculations for neu- 
tral defects, as this can lead to transition levels calculated in 
individual supercells (without scaling) to appear to lie inside 
the gap when they should lie outside and vice versa. Indeed, 
electrostatic errors are still present even for physically rea- 
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sonable cases, such as Vp and (probably) Vi n , due to defect 
states which enter the gap as the cell size grows. They are 
even present and may remain significant for defects such as 
the neutral Pi„ which has all of its defect levels within the gap 
in all supercells. This is because neutral defects in crystalline 
solids still have higher order charge multipole moments, es- 
pecially if the system is made up of more than one type of 
atom with different electronegativities. We have shown that 
for unrelaxed formation energies these errors scale as the in- 
verse cube of the supercell size (L~ 3 ). For relaxed formation 
energies this is almost certainly the case also, although we 
do find possible indications that the leading non-linear error 
term may sometimes scale as the inverse square (L~ 2 ). Fur- 
ther work with more defects is needed to confirm or defini- 
tively rule out this possibility. It could clearly be answered 
by removing the 8 atom cell from the scaling and replacing 
it with the 1000 atom cell, but that must wait for improved 
computing facilities. In the meantime, an answer may be ob- 
tained from scaling studies of formation energies for charged 
defects 12 since both the electrostatic errors themselves and the 
cross terms between them and the elastic errors will then be 
stronger. 

To summarize: the use of large (up to 500 or occasionally 
even 1000 atom) supercells with finite size scaling has been 
shown to be a very promising route around the errors which 
arise in the use of the supercell approximation to calculate for- 



mation energies of defects in III-V (and other) semiconduc- 
tors. Errors scale with a linear plus a higher order term, most 
probably cubic. We have also found several instances where 
scaling recovers physically relevant results that are even hid- 
den in calculations on the 512 atom cell: formation energies 
which are wrong by ~ A e V and defect levels which appear in- 
side the valence or conduction bands in supercell calculations 
when they should actually lie inside the band gap. 
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